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Abstract 



■ We report large systematic errors in Monte Carlo simulations of the tri- 

■ critical Blume-Capel model using single spin Metropolis updating. The error, 
manifest as a 20% asymmetry in the magnetisation distribution, is traced 

^ : 

0^ . to the interplay between strong triplet correlations in the shift register ran- 



dom number generator and the large tricritical clusters. The effect of these 



. correlations is visible only when the system volume is a multiple of the ran- 

O ■ dom number generator lag parameter. No such effects are observed in related 



models. 



I. INTRODUCTION 



Interest in random number generator (RNG) deficiencies has recently been aroused by 
the work of Ferrenberg et al |T|, who found that the widely employed "R250" generator 
produces systematic errors in Monte Carlo simulations using cluster updating algorithms. 
Their high precision measurements of the specific heat and the internal energy of the critical 
two dimensional Ising model differed from the exact values by up to 0.05% [T|]. 

These observations prompted a number of detailed studies of systematic errors in cluster 
algorithms and other "depth first" algorithms, such as those employed for simulating self 
avoiding walks |@|- 0. The consensus emerging from these studies was that pseudo random 
numbers (PRNs) generated by generalized feedback shift register (GFSR) algorithms 

of the type 

Xn = Xn-p®Xn-q (© = cxclusive-or) , (1) 

are generally "bad" when used in depth first algorithms, even though they pass many of 
the standard RNG tests and have very long periods [Q. However for local (single spin 
flip) Metropolis updating, good results were obtained using the GFSR generators. It was 
therefore concluded that this type of algorithm is relatively insensitive to correlations and 
thus safe to use. 

In this paper, we show that in certain circumstances large systematic errors can also arise 
when using GFSRs generators in conjunction with local Metropolis updating. We demon- 
strate this for the case of the tricritical Blume-Capel model and argue that the observed 
errors originate from certain extremely large triplet correlations in the random number se- 
quences. 

II. RESULTS 

There have been a number of recent suggestions in the literature that the problems 
observed with the GFSR generators may originate from triplet correlations [0^]. In view 
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of this we have made a study of three point correlations in the GFSR generators and have 
indeed found an extremely large effect in the triplet. For this triplet the 

three point average (X„ ■ ■ X„_g) takes the value 0.107 instead of (1/2)^ = 0.125. This 
effect is displayed in figure |1| where we plot the triplets (X„X„_fcX„_25o) as a function of k, 
for 1 < k < 249. The results were obtained from a sequence of 10® PRNs using the magic 
number pair (250, 147) corresponding to the R250 generator. In fact the same effect, with 
the same magnitude, is observable for other magic number pairs {p, q) recommended in the 
literature §, viz (521,168), (1279,216), (1279,418), (3217,67), (3217,576), (4423,1393), 
(4423,2098), (9689,84), (9689,471), (9689,1836), (9689,2444), (9689,4187). 

It is thus conceivable that triplet correlations are responsible for the systematic errors 
observed in "depth first" algorithms, which principally rely on long uncorrelated sequences of 
random numbers. However, as mentioned above, it has hitherto been generally believed that 
local Metropolis updating schemes are insensitive to correlations in the PRN sequences and 
give results at least better than linear congruential RNGs ||1],0. In fact we have found that 
in one system, the tricritical Blume-Capel model, triplet correlations can lead to extremely 
large (~ 20%) errors, even when using a Metropolis single site updating scheme. 

The Blume-Capel model is a spin 1 Ising model with the Hamiltonian 

nBC = -J2 S^S,+AJ2S^^ (2) 

<ij> i 

where spins Sj take the values = 0, ±1 and the sum < ij > runs over all nearest neighbor 
pairs. We have simulated this model on an FCC lattice for system size = (5 x 5 x 5 • 
4) = 500 spins. Spins were updated serially, each update requiring two random numbers: 
one to choose the new spin, and one for the Metropolis step. The random numbers were 
generated by the R250 RNG, having the magic number pair (p=250, q=147). Figure |^(a) 
shows the measured probability distribution of the magnetisation M = J^i^i at the near- 
tricritical parameters ksT = 3.0, A = 5.67 (circles). Notwithstanding the symmetry of 
the Hamiltonian under inversion of all spins {{si} the distribution is clearly 

asymmetric with respect to the "plus" and "minus" phases. We find, however, that if spins 
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are updated randomly instead of serially, the distribution becomes symmetric [figure 0(b)]. 

We believe that the observed asymmetry in P{M) is related to the fact that for system 
size V = 500, the larger generator lag p = 250 is a factor of the system volume, so that 
a given spin is always updated with a random number having the same relative position 
within successive sequences of 250 PRNs. This view is supported by the following evidence: 

• We find that the asymmetry in P{M) disappears when other system sizes are simu- 
lated, e.g. (14 X 14 X 14 ■ 4) = 1372. 

• The asymmetry disappears for the V = 500 system size when one additional PRN is 
drawn and discarded after each Monte Carlo step (i.e. after a whole sweep through 
the lattice). 

• The asymmetry disappears for the V = 500 system size when the R250 generator is 
substituted by the "R1279", a GFSR RNG with the magic number pair (1279,216). 

• The asymmetry reappears with the R1279, for the V = 500 system size when 279 
additional PRNs are drawn and discarded after each Monte Carlo step. 

We further believe that the asymmetry in P{M) is a manifestation of the aforemen- 
tioned triplet correlations in the PRN sequence (X„, X„_p, X„_q). To substantiate this we 
performed a simulation using the lagged Fibonacci RNG ^ 

Xn = Xn-p + Xn-g mod 2^^ with (p, q) = (1279, 216), (3) 

which is formally similar to the "R1279" , but does not display the triplet correlation. Even 
after discarding 279 PRNs after each cycle, P{M) remains symmetric with this generator. 

It is interesting to observe that a number of further factors influence the size of the 
asymmetry effect. One is the proximity to the tricritical point. If we choose A = 0, /c^T = 
6.8, corresponding to the critical point of the spin 1 Ising model, the effect weakens, but 
remains clearly discernible (not shown). We also note that the asymmetry effect is not 
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influenced by the dimensionality or the type of lattice: Simulations of the L = 10 x 10 Blume- 
Capel model on a square lattice (with 50 additional random numbers discarded between 
each sweep) exhibit the asymmetry with a similar magnitude to that observed in the three- 
dimensional case. 

The effect does, however, seem to depend on the symmetry of the Hamiltonian of the 
model studied. To demonstrate this we have also performed simulations on the three state 
Potts model near its critical point. This is a lattice spin model with the Hamiltonian 

Up = - ^s.,s,, (4) 

<ij> 

where 5ij denotes the Kronecker symbol and Sj takes again the values Sj = 0, ±1. The sole 
alteration required to the Blume-Capel program to simulate this model, is the calculation of 
the energy change in the Metropolis step. In figure ^ we show the magnetisation distribution 
P[M) for the three state Potts model on a = 500 FCC lattice, for the near-critical 
temperature fc^T = 3.9. No asymmetry between the "plus" and "minus" phases is evident. 
Similar simulations on the 2D simple Ising model also reveal no asymmetry. 

It is thus clear that the observed asymmetry in the tricritical Blume-Capel model is the 
result of a subtle interplay between the tricritical clusters and PRN triplet correlations. The 
fact that one of the three phases in the Blume-Capel model (the "0" phase) is distinguishable 
from the other two at the tricritical point seems to be crucial. From the joint distribution 
p(M, Q), with Q = Y.i^l (figure ^), one sees that there is no direct phase space path between 
the "plus" and "minus" phases. In order to migrate from one of these phases to the another, 
the system must pass through the "0" phase. It seems that the tricritical clusters of "0"s are 
very sensitive to correlations in the PRNs, and that "resonant enhancement" of the triplet 
correlations occurs when the system size is a multiple of the larger generator lag. The critical 
clusters in the Potts model are much less susceptible to these specific correlations. 
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III. CONCLUSIONS 



To summarize, we have demonstrated that the interplay between the physical properties 
of a system (principally the system size) and random number generator triplet correlations 
can lead to strong systematic errors in simulations even when simple Metropolis updating is 
used. Such "resonant enhancement" of PRN correlations can be very specific to a particular 
situation, as in our case of the tricritical Blume-Capel model, and not be observable at all 
in other, albeit very similar systems. Our observation therefore supports the view that a 
random number generator is only "good" in conjunction with the specific application for 
which it has been tested - every new application is also a new pseudo random number test, 
where known and unknown correlations may gain importance in an uncontrolled way. 
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and to P. Nielaba for practical help with some of the tests. Some of our simulations were 
carried out using externally provided versions of the R250, e.g. on the CRAY-YMP at the 
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FIGURES 

FIG. 1. The triplets (X„X„_fcX„_25o) as a function of 1 < A; < 249, generated from a sequence 
of 10® PRNs using the generator of equation |l] with p = 250, q = 147. 

FIG. 2. (a) Magnetisation distribution of the Blume Capel model on the V = 500 site FCC 
lattice at ksT = 3,A = 5.67, obtained using serial updating with the R250 generator, (b) The 
same distribution obtained using random site updating. 

FIG. 3. Magnetisation distribution of the three state Potts model on the V = 500 FCC lattice 
at ksT = 3.9, obtained with PRNs generated by the R250 RNG. 

FIG. 4. The joint distribution P{M, Q) of the near-tricritical Blume-Capel model. Also shown 
on the grid base is a contour plot of the histogram surface. 
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FIG. 1. The triplets {XnXn-k^n-25o) as a function of 1 < A; < 249, generated from a sequence 
of 10^ PRNs using the generator of equation 1 with p = 250, q = 147. 



8 




1.5 



1.0 



0.5 



0.0 



■1.0 



-0.5 



random updating 




0.0 

M 



0.5 



1.0 



FIG. 2. (a) Magnetisation distribution of the Blume Capel model on the V = 500 site FCC 
lattice at ksT = 3,A = 5.67, obtained using serial updating with the R250 generator, (b) The 
same distribution obtained using random site updating. 
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FIG. 4. The joint distribution P{M, Q) of the near-tricritical Blume-Capel model. Also shown 
on the grid base is a contour plot of the histogram surface. 
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